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ABSTRACT 

We present a theoretical model for primordial star formation. First we describe the structure of the 
initial gas cores as virialized, quasi-hydrostatic objects in accord with recent high resolution numerical 
^ studies. The accretion rate can then be related to characteristic densities and temperatures that are 

set by the cooling properties of molecular hydrogen. We allow for rotation of the gas core, assuming 
angular momentum conservation inside the sonic point of the flow. In the typical case, most mass then 
reaches the star via an accretion disk. The structure of the inner region of this disk is described with the 
standard theory of viscous disks, but with allowance for the substantial energies absorbed in ionizing and 
. dissociating the gas. The size of the protostar and its luminosity depend upon the accretion rate, the 

energetics of the accreting gas, and the ability of the radiation to escape from the stellar accretion shock. 
We combine these models for the infall rate, inner disk structure, and protostellar evolution to predict 
the radiation field that is the basis for radiative feedback processes acting against infall (Paper II). For 
realistic initial angular momenta, the photosphere of the protostar is much smaller and hotter than in 
the spherical case, leading to stronger radiative feedback at earlier stages in the evolution. In particular, 
once the star is older than its Kelvin-Hclmholtz time, contraction towards the main sequence causes a 
rapid increase in ionizing and far-ultraviolet luminosity at masses ~ 30M Q in the fiducial case. Since 
the cores out of which the first stars formed were much more massive than 30Mq and since feedback 
is dynamically unimportant at lower masses, we conclude that the first stars should have had masses 
> 3OM . 

Subject headings: cosmology: theory — early universe — galaxies: formation — stars: formation 

6 . 

The formation of the first stars marks the beginning of the long saga of galaxy formation and evolution. These stars, 
■ if relatively massive, are thought to bring the Universe out of the "dark ages" , heating and reionizing the intergalactic 
medium, and enriching their surroundings with metals. The stellar mass determines whether or not supernovae occur, 
and if so the energy of the explosion, composition of the ejecta, and type of remnant. Very energetic supernovae may 
produce gamma-ray bursts (GRBs). The first stars may have had a profound influence on the formation of supermassive 
rN \ black holes, globular clusters, and proto-galactic fragments. For all these reasons, as well as simple curiosity about the 
. very first stars in the Universe, we wish to understand the process of primordial star formation and how, in particular, 
the formation process may control the resulting stellar mass. 

When and how quickly the universe was reionized depends on the dumpiness of the intergalactic medium (IGM), the 
time of first star formation, its efficiency, and the ionizing luminosity per stellar baryon (e.g. Madau, Haardt, & Rees 
1999). The latter depends on the initial mass function (IMF) of the stellar population (Tumlinson & Shull 2000; Bromm, 
Kudritzki, & Loeb 2001; Ciardi et al. 2001; Schaerer 2002), although Bromm et al. (2001) showed that once the stars 
all have masses > 300M Q , the spectral luminosity per unit stellar mass becomes almost independent of the IMF. The 
hardness of the radiation field is also somewhat sensitive to the IMF (e.g. Tumlinson & Shull 2000; Schaerer 2002), so that 
He reionization can be affected. Significant reionization may also result from a population of supermassive black holes, 
making it more difficult to draw conclusions on the nature of the stellar population. The recent microwave background 
polarization results reported by Kogut et al. (2003) imply a relatively early epoch of reionization and a very top-heavy 
IMF for models invoking reionization due to stars. 

As with contemporary stars, the bolometric luminosities and the total energy output of the first stars depend on 
their mass. Thus, along with the overall efficiency of star formation, the IMF has implications for the contribution the 
primordial stars make to the near infra-red background (Santos, Bromm, & Kamionkowski 2002; Salvaterra & Ferrara 
2003). Unfortunately the observational determination of this background is hampered by the zodiacal foreground emission. 
Clusters of primordial stars may potentially be detected by the next generation of space-based near infra-red observatories, 
perhaps showing up in the number count versus flux distribution of very faint sources, or in the fluctuations of the 
background (Magliocchetti, Salvaterra, & Ferrara 2003). 

The metal enrichment of the IGM from core-collapse supernovae depends on the mass of the stellar progenitors. Massive 
primordial stars are thought to have much smaller mass-loss rates than their present-day cousins (Barrafe, Heger, & 
Woosley 2001) so the mass at core collapse may be quite similar to the initial mass. For supernova progenitor masses 
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> 260M Q , collapse is thought to proceed directly to a black hole, as the temperatures resulting from the pair-instability 
collapse are great enough to photodisintegrate nuclei, which counters explosive oxygen and silicon burning (Fryer, Woosley, 
& Heger 2001; Heger & Woosley 2002). There is expected to be relatively little mass ejection and metal enrichment from 
such supernovae. For smaller masses, down to about <~ 140M©, explosive O and Si burning is able to completely disrupt 
the star, leaving no remnant and ejecting large quantities of heavy elements. Progenitors from ~ 40 to 140Mq form 
black holes, with relatively inefficient metal ejection, while less massive supernovae form neutron stars, with more normal 
enrichment rates (however, see Umeda & Nomoto 2003). In principle, metallicity determinations from high rcdshift 
absorption line systems or from very metal poor local stars (e.g. Aoki et al. 2002; Christlieb et al. 2002) can constrain 
the IMF of the early generations of stars. 

A phenomenon that may be related to supernovae and hypernovae is the production of relativistic ejecta and thereafter 
a gamma-ray burst (Woosley 1993; Paczynski 1998; MacFadyen & Woosley 1999; Tan, Matzner & McKee 2001; Matzncr 
2003). The amount of relativistic ejecta accelerated by a blast wave in the stellar atmosphere is a sensitive function 
of the progenitor mass, density structure, and explosion energy (Tan et al. 2001). Massive stars that mostly collapse 
to a black hole and impart a large energy to the small fraction of ejected mass, are excellent accelerators of relativistic 
ejecta. However, for efficient gamma-ray production this ejecta should cither interact with a relatively dense circumstellar 
envelope, or be subject to internal shocks. If GRBs and associated afterglows are found at high redshift they will provide 
excellent probes of the IGM, for example through the study of 21 cm absorption lines (Furlanetto & Loeb 2002). 

The impact of the first stars on the heating and metal enrichment (and thus cooling) of the IGM may have important 
implications for the formation of other objects, such as supermassive black holes, globular clusters, and galaxies. If the 
first stars were sufficiently massive to collapse efficiently into black holes and their feedback on the surrounding gas was 
relatively weak, then they might become the seeds that grow to millions of solar masses and more. Other formation 
scenarios are possible (e.g. Rees 1984). One alternative requires direct collapse of gas to a black hole without a nuclear- 
burning stellar phase (e.g. Baumgarte & Shapiro 1999), and applies to "stars" with m» > 10 6 M Q (Zeldovich & Novikov 
1971). However, concentrating such a large mass of gas into a hydrostatic structure, without first forming less massive 
protostars, would be difficult for the very first generation of stars since molecular cooling acts to reduce the Jeans mass (e.g. 
Abel, Bryan, & Norman 2002, hereafter ABN). This scenario could perhaps be resuscitated for subsequent generations of 
star formation if the FUV background radiation produced by earlier generations of star formation were large enough to 
prevent H2 formation (Bromm & Loeb 2003). Other formation mechanisms have their origins in a very dense star cluster, 
composed either of nuclear-burning or compact stars. The ability of such clusters to form in the early universe is likely to 
be influenced by the nature of primordial stars; for example, if primordial stars are relatively massive, then the resulting 
feedback may inhibit further star formation in the vicinity (Omukai & Nishi 1999). The same considerations apply to 
early globular cluster formation. The dominant feedback on cosmological scales that affects galaxy formation is likely to 
be via the intensity of the FUV background. If this is very high, then by suppressing the H 2 abundance, the collapse of 
baryons into "mini-halos" would be inhibited. 

A description of primordial star formation requires relatively simple physics and chemistry, mainly because, by definition 
there are no complicated feedback processes from earlier stellar generations. For most early generations of stars, the only 
feedback variable will be the intensity of the FUV background, with perhaps some additional influence from hard X-rays 
(Cen 2003). In contrast to the present-day case, dynamically-important magnetic fields and dust grains are probably 
not present during the initial stages of collapse. The initial conditions for star formation are well-defined cosmological 
fluctuations. Given this relative simplicity, we have greater confidence in the applicability of the results of numerical 
simulations that have followed the collapse of cosmological scale perturbations down to almost stellar dimensions (ABN; 
Bromm, Coppi, & Larson 1999; 2002, hereafter BCL). This confidence is strengthened by the fact that it appears to be 
the microphysics of H 2 cooling that determines the types of baryonic structures that are formed, and not, for example, 
the details of the initial power spectrum of fluctuations in dark matter density. Molecule formation was predicted to occur 
in the presence of residual electrons left over from recombination via H~, (McDowell 1961; Hirasawa, Aizu, & Taketani 
1969). Including this physics, the results of the recent numerical simulations suggest that the initial gas cores out of 
which stars form are quite massive, M coro <~ 100 — 1000M Q , in agreement with earlier theoretical expectations (Yoneyama 
1972; Hutchins 1976; Silk 1977; Carlberg 1981). However, at very high densities three-body molecule formation becomes 
efficient, and it was speculated that the increased cooling would lead to much smaller Jeans masses <~ M (Palla, Salpeter, 
& Stahler 1983). The simulations indicate that this process does occur, but produces only a single initially sub-solar mass 
protostar at the center of the much larger core. Accretion then builds up the mass of the star (Omukai & Nishi 1998, 
Ripamonti et al. 2002). 

The goal of this paper is to model the growth of the protostar from very small masses to large, and to determine 
the total energy output from the star. Here we do not consider how feedback processes may inhibit the accretion to the 
protostar, deferring this important question to other papers. In §2 we briefly review the stages of the collapse immediately 
prior to protostellar core formation. We then describe the accretion rate to and the structure of the infall envelope around 
the protostellar core and disk (§§3, 4). We model the inner accretion disk in §5 and the evolution of the protostar in 
§6. The results for the bolometric and ionizing luminosities are presented in §7. Appendix A describes the accretion of 
energy by a protostar and its disk. General results on the structure and luminosity of an accreting protostar are given in 
Appendix B. Radiative feedback processes are considered in Paper II, while magnetic field generation and hydromagnetic 
winds have been modeled by Tan & Blackman (2004). 
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2. OVERVIEW OF BARYONIC STRUCTURE FORMATION: FROM COSMOLOGICAL FLUCTUATIONS TO STAR-FORMING CORES 

A simple analytic picture of baryonic structure formation is the following (e.g. Peebles 1993; Tegmark et al. 1997; 
Madau 2002). Recombination at a redshift z ~ 1200 heralds the start of the "dark ages" . Thermal equilibration of matter 
and radiation is maintained via the residual free electrons until z ~ 160. Until this point the cosmological Jeans mass (gas 
plus dark matter), Mj oc (T/p 1 / 3 ) 3 / 2 , is therefore independent of z and has typical value ~ 1O 5 M , similar to the scale 
of globular clusters. At lower redshifts, baryons cool adiabatically so that T oc (1 + z) 2 and Mj oc (1 + z) 3 / 2 . Gas collects 
in dark matter halos and cools further via H 2 line cooling. At some point the first luminous objects form, reionize the 
Universe to a temperature of <~ 10 4 K and raise the Jeans mass to galactic scales of <~ 10 9 ~ 10 [(1 + z rc _i on )/lO]~ 3 / 2 -M0. 

Numerical simulations (ABN, BCL) have largely confirmed this picture. Collapse is seen to proceed along filaments, 
with objects of characteristic mass Mj forming and rapidly merging to build-up somewhat larger masses. The simulations 
of ABN use adaptive mesh-refinement and follow cosmological perturbations in a relatively small volume from z — 100 to 
z ~ 20, where the first luminous object forms. These simulations are stopped when molecular lines, the main coolant at 
this stage, become optically thick. Spherically symmetric simulations (Omukai & Nishi 1998; Ripamonti et al. 2002) are 
able to follow the collapse all the way to the formation of a hydrostatic protostellar core. 

Summarizing the above numerical results: a pre-galactic halo (gas plus dark matter) of ~ 10 forms at the 

intersection of several filaments; a quasi-hydrostatic, gas-dominated core forms with mass ~ 4000M©, r ~ 10 pc, T <~ 
200 — 300 K (set by molecular cooling, with the mass fraction of molecules <~ 10~ 3 and their formation catalysed by free 
electrons via H~); gradual contraction of the core is driven by cooling in the dense interior; rapid 3-body H 2 formation 
occurs at u > 10 8 cm~ 3 so that most H is molecular by the time densities reach ~ 10 11 cm~ 3 , leading to strong cooling 
and supersonic inflow; and a hydrostatic core of mass ~ 5 x 10~ 3 M Q and radius ~ 14 R Q forms when the gas becomes 
optically thick to the cooling radiation (continuum cooling from collision-induced absorption). 

3. THE ACCRETION RATE: ISENTROPIC ACCRETION MODEL 

Consider a gas cloud that is initially in approximate hydrostatic equilibrium. We assume that this cloud undergoes 
gravitational collapse, forming one or more protostars. We focus on the gas that will go into a single protostar (or possibly 
a binary); we term this the protostellar core. We assume that the core is approximately spherical, with a baryonic density 
p(r) and mass M{r). Since the core is initially in approximate hydrostatic equilibrium, it will undergo an "inside-out" 
collapse (Shu 1977), in which the protostellar mass m* grows with time from a very small initial value to its final value. 

In their study of contemporary star formation, McKee & Tan (2002, 2003) assumed that the mass of any protostellar 
accretion disk was small compared to the stellar mass. The high accretion rates and possible absence of magnetic fields 
inferred for primordial star formation imply that, once most accretion is via a disk, then the rate limiting step for accretion 
is likely to be gravitational torques in the disk. Such torques become efficient only after the disk becomes relatively massive 
(disk mass m d ~ |m*) and self-gravitating (Adams, Ruden, & Shu 1989; Shu et al. 1990; Gammie 2001). Let 

77j*d = to, + m d = (1 + /d)m* (1) 

be the mass of the protostar and its associated accretion disk, where f d is the ratio of the disk to stellar mass. We adopt 
/d ~ 1/3 as a fiducial value. We allow for the possibility that only a fraction e*d of the mass in the core winds up in the 
protostar plus disk due to protostellar outflows. At any time after the collapse has commenced, we can then associate a 
core mass M{r) = m*d/c*d, and hence a density p(r), with the protostar and disk when its mass is m*<j- In contemporary 
star formation, the main driver for protostellar outflows is believed to be hydromagnetic winds (Konigl & Pudritz 2000; 
Shu et al. 2000). In this paper we assume that magnetic fields, and therefore protostellar outflows, are negligible in 
primordial star formation, so that e*d = 1; however, we shall keep e*<j in the equations for generality. (Tan & Blackman 
2004 consider the possibility of magnetic field generation by a disk dynamo, and the influence of the resulting outflow on 
e*d-) 

In general, the accretion rate to the protostar and any associated accretion disk can be expressed as 

A m * d iO\ 

m* d = 0* — - , (2) 

extending the expression given by Stahler, Shu, & Taam (1980) and McKee & Tan (2002, 2003) to the case in which 
the disk mass is a significant fraction of the protostellar mass. Here 0* is a numerical parameter of order unity and 
tfi — (37T/32G/3) 1 / 2 is the free-fall time measured at the location in the initial core where the interior baryonic mass is 
M = m,i,d/e*d- We assume that the mass fraction of dark matter is small, which is valid < 0.3 pc from the core center 
(ABN; note, however, that they assumed a baryon to dark matter ratio of 0.064, compared to 0.21 derived from the results 
reported by Spergel et al. 2003, so their results overestimate the effect of dark matter). On larger scales, the effect of 
dark matter is to reduce the free-fall time, or cquivalcntly increase the value of 0*, by factors of order unity relative to 
the purely baryonic case. 

Following McKee & Tan (2002), we now assume that the core out of which the star is forming is in a self-similar virial 
equilibrium with a density p(r) oc r p and a pressure P oc r~ kp . The pressure then obeys the polytropic relation 

P = K p p», (3) 

with K p = const and 7 P = kp/k p . In hydrostatic equilibrium, kp = 2(k p — 1), so that j p = 2(1 — l/k p ). If the gas changes 
in time, e.g., if it undergoes a reversible compression or rarefaction, then the change in pressure is related to that in the 
density by the adiabatic index 7 = 5\nP/S In p. If we define the parameter K by P = Kp 1 , then if is a function of 
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the entropy and is therefore termed the "entropy parameter" (McKee & Holliman 1999). If the adiabatic and polytropic 
indexes are equal (7 = 7 P ), then K = K p : a polytropic gas is isentropic and vice versa; furthermore, if the gas is initially 
polytropic, it will remain so under adiabatic compression or expansion. 

McLaughlin & Pudritz (1997) determined the accretion rate for a collapsing isentropic, singular polytropic sphere, and 
based on their results McKee & Tan (2002) showed that in this case 0* ~ 1.62 — 0A8k p . If the gas is not isentropic 
(7 7p); the value of 0* should not change significantly provided that the gas is convectively stable (-f p < 7) and is not 
stiff (7 < |). For the case of primordial star formation the assumption that the gas is isentropic should be quite good: 
BCL and ABN show that H 2 cooling causes the gas to accumulate at a density of order 10 4 cm' 3 and a temperature of a 
few hundred K, which sets the value of the entropy; Omukai & Nishi (1998) have shown that the subsequent compression 
of the gas obeys 7 ~ 1.09, leading to an isentropic gas distribution with 7 ~ 7 P ~ 1.1. 

McKee & Tan (2002, 2003) expressed the protostellar accretion rate in terms of the mass of the star and the pressure 
at the surface of the core out of which the star is forming. However, the accretion rate can equally well be expressed in 
terms of the entropy parameter K (Yahil 1983). From the equation of hydrostatic equilibrium, one can readily show that 
for an isentropic gas with 7 = j p , 



(3 - k p )k 3 P K 3 



4irG 3 M 2 



l/(4-3 7p ) 



(4) 



where M is the mass of the core interior to the point at which the density is p. Using equation (4) in equation (2), we 
then find 

(3-M4^ 3 i 1/[2(4 - 37p)1 
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where 
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(5) 
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For m^d oc M (actually, we are assuming m*d = M), this corresponds to a time dependence m*d oc fr 5 ^ 1- ^. Thus, 
whereas contemporary star formation accelerates (j p < 1 so that — j] > — McLaughlin & Pudritz 1997; McKee & 
Tan 2002, 2003), primordial star formation decelerates (4/3 > -y p > 1 so that j/[l - j] < 0) (Omukai & Nishi 1998). 

For 7 ~ 1.09 (Omukai & Nishi 1998), we have k p ~ 2.2, so that m*d oc M~ - 37 . The dependence of m*d on M is quite 
sensitive to 7„. Ripamonti et al. (2002) followed the growth of the protostar up to about 0.5M Q , and their result of 
m,j oc M -0-5 , i.e. j = —0.522, implies an effective value of j p = 1.114. We shall adopt an intermediate value 7 P = 1.1 
for the entire collapsing cloud. With this value we have k p = 20/9 and kp = 22/9. 

The normalization of the accretion rate then depends on two parameters, K and <j>*. First consider the entropy 
parameter K. As remarked above, ABN and BCL show that before the gas can collapse, it passes through a stage when 
its density is about 10 4 cm~ 3 , which is the density at which H2 is thcrmalized, and a temperature of about 200 K, 
which is the minimum temperature to which H2 can cool the gas. However, the core can also be supported by turbulent 
motions: in the simulation of ABN it is permeated by a cascade of weak shocks, whose velocity is about the sound speed 
(T. Abel, private communication, 2002). The total pressure is P = pc 2 = p{c 2 h + of urb ), where c is the effective sound 
speed, c t h is the thermal sound speed, and cr tur b is the 1-D velocity dispersion of the turbulent motions. In this case 
we take a — c t h/V3, so that the total pressure is a factor 4/3 higher than the thermal pressure. Equivalently we may 
define an effective temperature T e g = P/(nk), that is a factor 4/3 greater in this case. We consider primordial gas with 
X = 0.76 and Y = 0.24 so that in the atomic phase p = 1.22toh and Uho = 0.079nH- The isothermal sound speed is then 
cth = \JkT j p = 1.425(T/300K) 1 / 2 km s _1 . Normalizing the entropy parameter to these values gives 
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The parameter 0* depends on the type of collapse solution, which depends on the initial condition: Collapse from a 
singular isothermal sphere leads to the Shu (1977) solution for inside-out collapse with 0* = 0.663, whereas analogous 
collapse from a singular polytropic sphere with 7 P = 1.1 gives 0* ~ 0.55 (McKee & Tan 2003). The corresponding value 
of the accretion rate is substantially less than that which occurs in the collapse of a uniform core (Larson 1969; Penston 
1969); for the isothermal case, it is 48 times less (Hunter 1977). The simulations of ABN suggest that primordial star 
formation proceeds in a manner intermediate between these two extremes, but closer to the inside-out collapse: at the 
time of protostar formation the envelope is contracting at about a third of the sound speed. Hunter (1977) has described 
the family of isothermal solutions that span the range between the Shu solution and the Larson-Penston solution. In 
particular his model (lib) has collapse velocities of about v m ~ c t h/3 for the gas at the time of protostar formation. 
Compared to the Shu-type solution (eq. 5), the accretion rate is about a factor 2.6 larger. Assuming the same increase 
applies to the 7 P ~ 1.1 case, we find (from eq. 5) 

-3/7 



m* d = 0.026e* d K' 15 / 7 (JL 



Mr. 



(9) 
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(10) 

>3/4 



(1 + .U) 3/J 

Note that the stellar accretion rate and the mass flow rate through the inner accretion disk is a factor (1 + fd)~ l 
smaller than this. The accretion rate declines in time as i~ 3//10 . The time required to build up a given stellar mass is 

/ \ 10 /7 

«. = 27(1 + / d rvrv- i5 / 7 ^ j y , (n) 

The lifetimes, t\a e , of very massive stars have been estimated to be ~ 2 Myr (Schaerer 2002). The condition t* < t\a e 
implies that a maximum stellar mass of about 2000M Q could be accumulated in the absence of any feedback processes. 
On these scales the effect of dark matter becomes important, so that the formation timescale for a given baryonic mass 
becomes shorter, raising the above mass limit by a modest factor. 

The accretion rate for this model is shown in Figure 1, along with the rates predicted by Ripamonti et al. (2002) and 
Omukai & Nishi (1998). Omukai & Nishi assumed that the collapse would follow the Larson-Penston solution, which is 
appropriate for an initially uniform cloud. In the Larson-Penston solution, the infall velocity is about three times the 
sound speed at the point of protostellar core formation. However, the calculations of Abel et al. (2000, 2002) show that 
the protostellar core is highly centrally concentrated and has only subsonic infall motions. This is to be expected, since 
it is evolving quasistatically due to the slow radiation losses from the trace amounts of molecular hydrogen present in 
the primordial gas. In order to fit their observed accretion rates with an analytic model, ON were forced to adopt an 
artificially small value of the entropy parameter K, because of their use of the Larson-Penston solution. Equation (10) 
is a more accurate description of how the normalization of the accretion rate depends on the properties of the pre-stellar 
core. 
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Fig. 1. — Mass accretion rate onto the protostar and disk as a function of the current collapsed mass, m,^ (we assume e t( j = 1, so that 
m t£ j = M). Solid line fiducial model (with K' = 1) from eq. (10). Dotted line from Omukai & Nishi (1998). Dashed line is analytic result 
from Ripamonti et al. (2002). Long-dashed line is the constant accretion rate used in protostellar evolution models of Stahler et al. (1986) 
and Omukai & Palla (2001). Dot-dashed line is the settling inflow rate at the final stage of the simulation of ABN, now as a function of the 
enclosed mass. Note that the decline of this rate at small masses is due to the lack of the full set of high density cooling processes in the 
simulation of ABN. 



3.1. The Mass Scale of Primordial Star Formation 
The mass of the protostellar cloud can be inferred from equation (4): 

10 4 cm 



M = 543K' 372 



n H 



7/20 



M G 



(12) 



where M is the mass of gas at densities greater than tih under the assumption that most of the mass is polytropic. This 
mass is somewhat less than the maximum possible mass that can be gravitationally stable for these values of 7 P and tih, 
as is generally true for polytropes (see Fig. 6 in McKee & Holliman 1999). 
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We can estimate the mass of the molecular core in this cloud by equating the formation time for H2 to the dynamical 
time, which we take to be t*. Rapid H2 formation can occur only at high densities via a three-body process. The rate 
constant, fc 3b , for H + H + H -> H 2 + H used in the simulation of ABN was k 3h = 1.3 x lCT 32 (T/300 K)~ - 38 cm 6 s" 1 for 
T < 300 K (Orel 1987) and fc 3b = 1.3 x 10~ 32 (T/300 K) _1 cm 6 s" 1 for T > 300 K (Palla et al. 1983). The characteristic 
time for the gas to become molecular is then 

*h 2 = ~^r~~ — 2" = 1-22 x l^K' ^ y , (13) 

Equating this to the star formation time t* (which may be inferred from equations 11 and 12 with to* = e*<jM/(l + fd)), 
we find that the gas becomes molecular at a density nn = 2.8 x 10 n K' 5 / 7 cm -3 , corresponding to a mass m. t d = 
l.34e„ d K /5/4 M Q . Omukai & Nishi's (1998) and ABN's results show that the gas becomes molecular at T ~ 600 K, so 
it has a slightly different entropy parameter (K' ~ 0.4) than the gas at tih = 10 4 cm -3 . We conclude that the central 
0.4M Q of the initial gas cloud is molecular, which is close to the results of Omukain & Nishi (1998) and ABN (as shown 
in the final timcstcp of their simulation). 



4. ROTATING INFALL 

The collapsing clouds formed in the simulations of ABN and BCL are rotating. When the initial angular momentum 
is very high, then the structure may form a rotationally-supported disk on quite large scales ~ 10 3 M Q (Bromm et al. 
1999). It then fragments into protostellar cores with masses ~ IOOMq. For more typical cosmological initial conditions, 
ABN find that the first "molecular cloud" that forms in their simulation does not fragment. It is rotating, but never (on 
the scales resolved ~ 20 AU) fast enough to be completely supported by rotation. Their figure 4b shows that the radial 
mass- weighted rotational velocities are about half of Keplerian for a broad range of the enclosed mass. These rotational 
velocities are therefore also about the same magnitude as the sound speed. ABN attribute the angular momentum 
transport needed to maintain rotational velocities v Iot ~ |^Kep to weak shocks in the flow. They were unable to continue 
their calculation to the point that supersonic infall commenced. Once such infall starts, it will be difficult to transfer 
angular momentum from inside the sonic point to outside (indeed, in the absence of shocks in the flow, it would be 
impossible). We idealize the angular momentum transport in the infalling gas by assuming that u ro t <x «Kcp outside 
the sonic point (specific angular momentum j oc r 1 / 2 ) and that the angular momentum is constant inside. The angular 
momentum of the gas that accretes onto the star-disk system is then characterized by the parameter 



, _ ^rot(r sp ) Vrot{r S p, , . 

/Kep " v Kcp (r sp ) (GM sp /r sp ) 1/2' 



where r sp is the radius of the sonic point and M sp is the mass interior to it. Averaging over spherical shells, ABN found 
w rot ~ 0.5wko P j independent of r, so we adopt 0.5 as a fiducial value for /xep- 

Conservation of angular momentum from the sonic point to the outer radius of the accretion disk, rd, implies 



where the mass interior to r d is to^ = (1 + /<j)m*. We shall typically assume that no mass is diverted from the inflow 
between r sp and r d , so that in the above equation, M sp = m^d- In order to evaluate M sp (r sp ), we note that the mass-radius 
relation of the equilibrium core is 



M 



47T 



1-7 P 



3-kJ ( G ) 



r 4-3 7p 



1/(2-7*0 



^ 980 f-^J K' W / 9 M Q (16) 

from equation (4) together with the relation M — Anr 3 p/ (3 — k p ). The mass inside the sonic point is made up of the core 
of mass TOoc 3 h t/G (with to = 2.577 for the isothermal settling solution vs. to = 0.975 for the Shu solution) and the infall 
envelope, which has a mass ~ c^ h t/G for both solutions (see Fig. 3 in Hunter 1977; most of the mass is near ( ~ 1). Thus 
in the isothermal case the mass within the sonic point radius of the settling solution is about (2.577+ l)/(0.975 + 1) ~ 1.8 
times greater than in the equilibrium state. We assume a similar increase applies to the 7 P ~ 1.1 case. If the protostar 
generates a bipolar outflow that sweeps up material, then the mass inside the sonic point will be smaller by a factor of 
about e*d- Using this as the efficiency at this stage, and under the assumption that no mass is diverted from the inflow 
between r sp and r<j, we then have to*^ ~ M sp = 1.8e*dM(r sp ). Numerical evaluation of equation (15) then yields 

rd 3 - 44 (t?) 2 ^ 77 (S) 9/7 ^ 10/7 AU - 

Equation (17) defines the outer radius of the disk. Material falls on to this disk at all radii r* < r < r d (Figure 2). 
There is also some direct accretion to the star of material that had very little initial angular momentum. We follow Ulrich 
(1976) in describing the density distribution of this freely-falling and rotating accretion envelope. Note that Ulrich and 
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those who followed (Cassen & Moosman 1981; Terebey, Shu, & Cassen 1984) assumed that the rotation was initially 
uniform, so applying this solution to a turbulent, differentially rotating cloud is necessarily approximate. The fact that 
the characteristic disk scale is generally much greater than the stellar radius (which is always less than a few hundred 
solar radii, §6) implies that most accretion proceeds via a disk. Provided that > r*, the mass accreting directly onto 
the disk is 



(Adams & Shu 1986). Once the disk is large compared to the star, direct accretion onto the star accounts for only a 
fraction r*/2rd of the total. 

The properties of protostellar disks supplied by inflow with the geometry shown in Figure 2 have been considered by 
Cassen & Moosman (1981). The impact of gas streamlines with the disk provides some viscosity to enable accretion 
through the disk, although this is not important for r <C r^. However, if the disk becomes sufficiently massive we expect 
the gravitational turbulent viscosity of clumps in a massive disk (Gammie 2001) or large scale (ra=l mode) instabilities 
(Adams et al. 1989; Shu ct al. 1990) to be the most important processes for driving inflow. Dynamo amplification of seed 
magnetic field (Tan & Blackman 2004) may lead to the development of magneto-rotational instability (MRI) (Balbus & 
Hawley 1991), which can provide a source of viscosity. Fragmentation of a gravitationally unstable disk may occur if the 
local thermal time of the disk becomes less than about half the orbital time (Gammie 2001), but this does not seem to 
occur in these disks (Tan & Blackman 2004). The structure of the inner accretion disk is considered in §5. 

One consequence of the geometry of rotating infall for initial protostellar core formation is that the optically thick 
stages of the collapse (with respect to both molecular line and continuum cooling, see §2) are reached at somewhat higher 
masses than in the case of spherical accretion. Since material can dissipate some of its energy in the disk, it has a smaller 
bulk kinetic energy when it reaches the star. It also has a higher temperature because of the dissipation of this energy, 
which has implications for the temperature of the post-accretion shock relaxation region and the protostellar evolution 
(see below). By changing the size and temperature of the photosphere around the star, disk accretion tends to create a 
hotter radiation field, which has an important effect on the feedback processes that occur (§7; Paper II). The density and 
a simple model of the optical depth near the star are shown in figure 3. 




(18) 



whereas that accreting directly onto the star is 




(19) 



2 



/ 



/ 



1.5 - 



/ / 
/ / 
/ / 
/ / 




0.5 1 1.5 











2 



Fig. 2. — Streamlines of the rotating, infalling envelope, projected onto the meridional plane. The protostar is at (0,0). The regions between 
each streamline and/or the axes, when rotated about the 2-axis, each contain 10% of the total accretion flow from this hemisphere. 
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Fig. 3. — Top panel: Density of flow relative to that at along radii at angles {dotted), ir/3 (solid), and 9n/20 (dashed) from the rotation 
axis. Note that the absolute densities at of these three cases are in the ratios 0.33:0.59:1.43 relative to the spherical case. Bottom panel: 
Characteristic optical depth (t') in traversing length scale r (assuming constant opacity), relative to that at r<j, for the same cases as above. 
Although ft is a more complicated function of p and T, this figure illustrates that the optical depth near the star is much reduced with respect 
to the spherical case and that under certain conditions when the envelope is already optically thick, it may become optically thin from the 
inside out, as the overall mass accretion rate declines. 



5. PROPERTIES OF THE INNER ACCRETION DISK 

Since much of the accretion onto the star goes through an accretion disk, it is necessary to determine the properties of 
this disk, particularly near the star where the disk emission can have significant feedback effects. Our objective in this 
section is to estimate the luminosity and approximate spectrum emitted from the inner region of the accretion disk. The 
luminosity from the inner part of the accretion disk can be comparable to that produced at the boundary layer. The 
structure of the accretion disk, and in particular its scale-height, is also of interest because it affects the size of the region 
where the accretion luminosity from the boundary layer is released, and thus the spectrum. The disk may be the site 
of dynamo action that can amplify weak seed magnetic fields and perhaps generate a hydromagnetic bipolar outflow, a 
common occurrence in contemporary star formation. For the primordial case this question has been addressed by Tan & 
Blackman (2004), where it is concluded that strong large-scale magnetic fields are likely to be produced if turbulence in 
the disk generates hclicity, particularly by the time the protostar is massive and contracting to the main sequence. Even 
in the absence of magnetic fields the disk may generate a radiatively driven outflow (Paper II), which depends somewhat 
on the structure of the disk. Finally the disk may act to shadow gas in the equatorial plane from the brunt of the stellar 
radiative feedback, thus influencing the efficiency of star formation. 

In order for the star to continue accreting, it must keep its rotation rate slower than the break-up rate. This is most 
likely achieved either by the action of magnetic fields coupled to the larger scale disk or an outflow, or by the excitation of 
density waves in the disk by the star (if it is rotating so quickly that it starts to assume a non-axisymmetric structure). We 
have assumed that these processes occur; furthermore, we assume that they are sufficiently effective that they maintain 
the overall rotation of the star to be much less than break-up, so that the rotational energy is negligible. 

To calculate the radial structure of the inner accretion disk at any given point in the evolution of the protostar, we 
assume that it is fed smoothly at a rate given by equation (10) and use the standard theory of steady, thin, viscous 
accretion disks, with a spatially constant viscosity parameter, a ss (Shakura & Sunyaev 1973; Frank, King, & Raine 1995). 
The viscosity is assumed to be a function of the total pressure, which is the sum of gas and radiation pressures. We have 
evaluated two cases with a ss = 0.01, 0.3 but regard the a ss = 0.01 case as the most appropriate for the very inner regions 
that are relevant here (see Tan & Blackman 2004 for a discussion on these choices and for the full results from both cases). 
This is a typical value for viscosity generated by the magneto-rotational instability (MRI) (Balbus & Hawley 1991, 1998), 
although the dispersion in quoted results from numerical studies is about two orders of magnitude. It must be stressed 
that the use of a spatially constant "a ss " viscosity model is motivated by its theoretical simplicity and convenience, and 
results of this modeling should be viewed as being only an approximate guide to reality. 
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The inner scale of the disk is set by r*, which we estimate in §6 below. We ignore the effects of energy injection from 
the star, deferring this to paper II. We again use the opacities of zero metallicity gas of Rogers & Iglesias (1992) and 
Iglesias & Rogers (1996) for T > 6000 K and of Lenzuni et al. (1991) for T < 6000 K (however, here the interest is only 
in the higher temperature regime). The method of solution assumes that the mean disk temperature, T, is much greater 
than the surface temperature, T c g, which is a good assumption for the inner regions that we consider, except in the early 
stages. 

The high accretion rates and large sizes of primordial protostars cause the ionization energy to have an important effect 
on the disk structure, particularly during the earlier stages of the evolution. The energy equation for the disk is worked 
out in Appendix A.l: 



87rr 3 



1 



(t) 



1/2 



777* d 

Anr dr \ 3 Cth 



<7 



(20) 



where e tn and ej are the thermal energy and dissociation/ionization energy per unit mass, averaged over the thickness 
of the disk. The first term on the right hand side of eq. (20) is the viscous dissipation per unit disk face area. Note 
that usually de^/dr < and dei/dr < 0. The effect of the thermal energy term is small for thin disks such as those 
we consider, so we drop it from our numerical calculations. However, the dissociation/ionization term, which is often 
neglected, can be an order of magnitude larger than thermal term and can have a major impact on the disk structure. 
We therefore retain this term. 

In this paper, we give an approximate solution for the disk based on the treatment of Frank et al. (1995). They 
integrate the radiative diffusion equation to obtain F ~ 4crT c 4 d /3r, where T C]d is the midplane temperature and r is the 
optical depth from the midplane to the surface. To solve the structure of the disk we divide it into a large number of 
radial zones, and start at the outer edge with an assumed ionization state that is almost completely neutral. We then 
solve the modified disk structure equations to find the temperature, density, ionization state of H and He (from the Saha 
equation), etc., and then advance the solution inwards given the new ionization state. Figure 4 shows three examples of 
disk structure for a ss = 0.01, with m* = 1, 10, 1OOM . At these masses the stellar sizes are taken to be r* = 100, 300, 4i? 
and the accretion rates in the inner disk — i.e., onto the star — are 777* = (17,6.4,2.4) x 1O _3 M yr _1 , respectively (eq. 
10). 

As we shall see below, at early times the large accretion rates in primordial protostars lead to large rates of energy 
absorption by ionization and dissociation (eq. 25), and this can regulate the midplane temperature to about 10 K in the 
early stages. The effect of including and not including the ionization energy is shown in the temperature plots of Figure 4. 

Our disk model breaks down when the predicted surface temperature is very cool (T < 5000 K), which occurs for low 
stellar masses. This break down occurs because in our thin disk model we approximate the opacity by its value at the 
midplane, whereas in fact the opacity increases rapidly with temperature for T < 10 4 K. However, since we are most 
interested in the radiative feedback from the star and disk at higher masses, we do not attempt to correct this limitation. 



6. PROTOSTELLAR EVOLUTION TO THE MAIN SEQUENCE 

The basic process we wish to model is the evolution of the protostar's radius as it grows in mass by accretion of material 
from an accretion disk or directly from the infalling envelope. Stahler et al. (1986) and Omukai & Palla (2001) considered 
a detailed model for the evolution of a protostar with spherically symmetric accretion and a constant accretion rate of 
4.4 x 1O~ 3 M yr _1 . Like these authors, we also define the protostellar radius to be the location of the accretion shock 
(see Appendix B). This size is important for determining the accretion luminosity of the protostar and it sets a lower 
limit for the photospheric radius of the protostar, which controls the spectrum of emitted radiation. 

In order to be able to follow the evolution of a protostar for an arbitrary accretion rate, we use the analytic approach 
developed by Nakano et al. (1995; 2000). They estimated the protostellar radius by equating the rate of accretion of 
energy by the protostar (which is worked out in Appendix A) with the rate of change of the energy of the protostar, 

2 



E 



+ m»e/. 



(21) 



The first term in this equation represents the sum of the thermal energy of the gas, the energy of the radiation and the 
gravitational energy, and it follows directly from the virial theorem. Here (3 is the mean ratio of gas pressure to total 
pressure and a g Gml/r* is the gravitational energy. We approximate the structure of the protostar as a polytrope of index 
77, so that a g = 3/(5 — 77) with 77 < 5. The second term in equation (21) is the mean ionization energy of the material in 
the star. We shall assume that the star is approximately fully ionized (ej ~ e/ TO = 16.8 cV/mn), which means that we 
cannot treat the earliest stages of protostellar evolution with our method. Differentiating equation (21) and equating it 
to the rate of energy accretion accretion given in equation (A10) evaluated just inside the postshock relaxation layer of 
the accretion shock (denoted r 2 ) yields 



dlnr* rfln/3 4 
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where /i 2 is the mean mass per particle at r 2 and E nuc is the rate at which nuclear reactions release energy in the protostar. 
This result is identical to that of Nakano et al. (2000), except that (1) equation (22) has the enthalpy 5fcT/2/j instead 
of the thermal energy 3fcT/2/x, (2) we allow the accretion flow to be optically thick, so that it is possible for T 2 to be 
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Fig. 4. — Protostellar disk structure for models with a as = 0.01 and m* = 1,10, 1OOM0, for which r* = 100, 300, 4/Jq (see §6) and 
rh„ = (17, 6.4, 2.4) X 10 — 3 Mq yr — 1 (see §3), respectively. From top to bottom the panels show (1) surface density, E; (2) ratio of scale-height 
to radius, h/r, and ratio of gas pressure to total pressure, f}\ (3) ionization fractions of H, He + , He 2+ ; and (4) central disk temperature, T Cj d 
(the dotted lines show results for when the ionization energy is neglected). 



much larger than the photospheric temperature, and (3) we explicitly allow for both direct accretion onto the star and 
for accretion via a disk. 

We now consider the evaluation of each of the terms in this equation. First, the term <iln/3/cilnm» becomes significant 
only when radiation pressure becomes important, which occurs only at relatively high masses. The appropriate polytropic 
index in this regime is n = 3, and in this case the mean value of [3 is the same as the central one. We calculate f3 assuming 
a standard Eddington stellar structure model with n = 3 applies to the protostar. 

Next consider the temperature behind the accretion shock, which is discussed in Appendix B. For accretion directly 
onto the star, the value of T2 in the optically thin case is given by equation (B5). As remarked in Appnedix B.l, we allow 
for the non-blackbody nature of the radiation field there by using the Saha equation to evaluate e/2 at T c g 1 2 ■ When the 
accretion is optically thick, we integrate through the radiative precursor to find T 2 , as described in Appendix B.2. In the 
optically thick case, e/2 is evaluated at T 2 . 

Disk accretion presents a more complicated problem. Gas reaches the surface of the star at a temperature T that we have 
estimated from standard thin disk theory (§5). However, as this gas joins the star, it will spread over the stellar surface, 
where it can lose heat. Integrating the radiative transfer equation through the disk shows that the mean temperature is 
related to the effective temperature by T oc t 1 / 4 ! 1 ^ (Shakura & Sunyaev 1973). For a constant opacity per unit mass, 
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r oc 1/A, where A is the area normal to the direction of propagation of the radiation; similarly, T c g oc 1/A 1 / 4 , so that 
altogether T oc 1/A 1 / 2 . For the disk, we identify the A as the area of the star covered by the disk, 47rrJcos^, where 
9d is the angle from the axis at which the disk intersects the star. We therefore adopt the approximation that for disk 
accretion, T 2 = T(r* ) (cos Od) 1 / 2 - In this paper, we estimate the height of the disk as 1.5 times the scaleheight; a more 
accurate estimate will be given in Paper II. 

To evaluate the luminosity L 2 , which is effectively the luminosity that can be transported radiatively, we follow Nakano 
et al. (1995) and approximate it as being equal to the value of the ZAMS luminosity (Schaerer 2002) of a star of the 
same mass. However, this is not a very good approximation at lower stellar masses, to* < IOMq, when, because of its 
expanded size, the protostar's internal luminosity is much smaller than the ZAMS value (Stahler ct al. 1986). In this 
regime we evaluate L 2 from analytic fits to the results of Omukai & Palla (2003, fig. 5), so that L 2 — for to* < 7M Q , 
then rises as L 2 /Lq ~ 390(to*/M q ) 2 — 2700to*/M q , which is used until L 2 — £zams- The idea behind this approach 
of specifying L 2 (m^,) is that the amount of radiative luminosity that a star can carry is determined by its structure. In 
a main sequence star, the rate of energy generation adjusts to give this luminosity; in a rapidly accreting protostar, the 
rate of energy generation, both nuclear and gravitational, can differ substantially from that which can be transported, 
and this difference can cause the star to cither swell or shrink. The evolution during the early stages is not particularly 
sensitive to this approximation because the accretion luminosity is much greater than the internal luminosity. 

The rate of nuclear energy generation, E nuc , is determined by both deuterium burning and the nuclear reactions that 
occur on the main sequence. We do not explicitly include the latter, but note that they prevent r* from shrinking below 
its main sequence value. Deuterium burning sets in when the central temperature T Ci * ~ 10 6 K. The central temperature 
is 

T c „=f3 c a T ^^ (23) 

K 7** 

where /x is the mean particle weight in the ionized interior and ax is a quantity of order unity (e.g. ay = 0.54,0.84 for 
n = 1.5,3, Chandrasekhar 1939). For m* < 10M Q , which is the regime where D-burning can be relatively important, 
/3 gas and (3 C are close to unity. The reaction rate rises sharply with temperature so that, while there is fuel available, this 
process can act as a thermostat and maintain an almost constant T Cj * (Stahler 1988). If the protostar reaches this central 
temperature at relatively low masses (as in contemporary star formation), then the luminosity produced by D burning 
can be large compared to that which can be transported radiatively through the star. The star becomes convective, which 
allows freshly accreted deuterium to be brought to the center for burning (we take D/H = 2.2 x 10 -5 ; Pettini & Bowen 
2001). The rate of D-burning can be calculated in this case using the model of Nakano et al. (2000) (see also McKee 
& Tan 2003). However, for realistic accretion rates (i.e. from cores with K' <~ 1), the protostar reaches the D-burning 
temperature only at relatively high masses (to* ~ 10M Q ) so that the star remains radiatively stable, D is quickly depleted 
from the stellar center, and the protostar can continue to evolve to hotter central temperatures, little affected by this 
nuclear energy generation. In this case we set E nuc = 10 4 Lq, independent of stellar mass, which approximates the more 
detailed calculations of Omukai & Palla (2003). This choice has only a minor effect on the stellar radius for most of 
the D-burning regime since usually E DUC <C L 2 , the internal protostellar luminosity (Fig. 6). Eventually the central 
temperature becomes hot enough to allow support via hydrogen burning. Some H-burning reactions start generating 
significant luminosity when T Cj * > 2 x 10 7 K (we set E nuc = 10 5 L & at this point, Omukai & Palla 2003). However, full 
support of the star only occurs once T Cj * ~ 10 8 K, which is hotter than in the contemporary case, leading to smaller 
main sequence radii (e.g. Schaerer 2002). In this way the protostar settles on the zero age main sequence, where it may 
continue to grow in mass if the accretion is on-going. 

Finally, the evolution of the polytropic index of the star needs to be accounted for. We take n to be constant during 
specific phases of the evolution, but allow for certain transitions as follows. By comparison to the results of Stahler et 
al. (1986), we find that an initial value of a g = 1.1 (corresponding to n = 2.3) is a reasonable description of the initial 
structure, while the star is younger than its Kelvin-Hclmholz time. If the central temperature becomes hot enough for 
D-burning and -E nuc ^ L 2 , so that the star becomes convective, then n = 1.5 (a g — 0.86) (this case is not realized in our 
models with K' = 1). In this case, the star eventually becomes radiative (n = 3 and a g = 1.5). The star also relaxes to 
this structure in the absence of convective D-burning, after the age of the star becomes greater than the current Kelvin- 
Helmholz time (this is the case realized in all the models presented in this paper). The expansion of the protostar after 
a Kelvin-Hclmholz time reflects the fact that it is relaxing to a globally more compact state, which releases gravitational 
energy and causes the outer layers of the star to expand. From comparison with the results of Omukai and Palla (2001) 
we set this expansion factor to be three. For protostars that evolve from a convective D-core burning state to a radiative 
state including D-shell burning, there is an analogous expansion of the outer layers, by about a factor of two (Palla & 
Stahler 1991). 

We take our initial condition to be a protostar of mass to* = 0.3Mq and radius r* = 30 Rq, which is a reasonable 
extrapolation of the numerical results of Ripamonti et al. (2002), who formed a core with a mass of O.O4M and a size 
of 10 12 cm ~ 14 Rq at the end of the core formation phase of their simulation (see also Omukai & Nishi 1998). In any 
case, since the gas in front of the accretion shock is optically thick at early times, the subsequent evolution is insensitive 
to the initial condition. 

A summary of the basic features of the protostellar evolution is as follows. The mean density and central temperature 
{T c ,* ~ 10 5 K) of the initial state are relatively low compared to typical stellar values. However, as the mass increases via 
accretion, self-gravity compresses the star and raises the central temperature. The rate of this contraction is limited by 
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the Kclvin-Hclmholz time; indeed rapid accretion of material on shorter timescales swells the star. The star is typically 
able to reach a mass ~ 10M Q by the point that its age is about equal to its instantaneous Kelvin-Hclmholz time. It 
also reaches the D-burning temperature at about this time: too high a mass for D-burning to create a convective core. 
Structural rearrangement to a radiative (n = 3) core temporarily expands the outer layers of the star by factors of a few. 
However, the rapidly increasing internal luminosity causes fast contraction (sec cq. 22) that is only halted once the star 
reaches the zero age main sequence, where the central temperature (T Cj * ~ 10 8 K) is hot enough for H-burning reactions to 
support the star. This occurs at stellar masses of about 100M Q . In our numerical implementation of the above evolution 
we allow for accretion from both a disk and directly to a star, the relative proportions of which depend on /kep, the size 
of the star, and the mass that has collapsed. 

Observe that because of the high accretion rates the total energy absorbed by dissociation and ionization (and the 
corresponding reduction in the protostar's luminosity that escapes to infinity) can be quite large: 



This process helps to keep the radiated luminosity sub-Eddington during the early stages of protostcllar evolution. In the 
later evolution of models of protostars accreting from realistic rotating cores, then equation (25) applies to the reduction 
in the accretion disk luminosity, as discussed in §5. 

A comparison of our simple semi-analytic model with the spherical accretion test case of Stahler et al. (1986) and 
Omukai & Palla (2001) is shown in Figure 5a. As discussed above, the expansion factor of the luminosity wave has been 
adjusted to give agreement with the models. We then applied the same model to the accretion rate of equation (10), but 
still with spherical symmetry. We show two examples of rotation for the isentropic model: the fiducial /k cp = 0.5 case 
and one with ,/kcp ten times smaller. 

Rotation causes the infall envelope to become optically thin at much smaller masses. This causes the photosphcric 
radius to be much closer to the stellar surface, leading to a hotter radiation field (§7). Also fk is reduced, because most 
material is processed through an accretion disk and so there is less bulk kinetic energy to advect, even if conditions are 
optically thick. However, the conditions at point "2", just behind the accretion shock are hotter in the disk case, so that 
overall more energy is advected into the star and its size is bigger. 



We are now in a position to calculate the total bolomctric luminosity of the protostar. This luminosity is the sum of the 
radiation from the stellar surface, the boundary layer, and that from the inner accretion disk; we somewhat arbitrarily 
include emission out to a distance lOr* in the latter. Emission from further out in the disk is not significant for the 
feedback processes we consider in Paper II. The luminosity from the boundary layer is comparable to, but almost always 
somewhat greater than, the luminosity from the inner disk. At high stellar masses the internal luminosity of the star 
becomes dominant. Figure 6a shows the evolution of the total luminosity and its various sub-components for the fiducial 
/Kep = 0.5, K' = 1, a ss = 0.01 case. The luminosites of the models with different rotational parameters are quite similar 
(Fig. 6b). The main differences are due to the different protostellar sizes (and thus accretion luminosities) and to the fact 
that our definition of the total luminosity only includes energy generated in the immediate vicinity of the star (r < 10r*) 
and not, for example, the outer accretion disk. This last point causes spherically-accreting protostars to have a higher 
total luminosity than disk-accreting ones, and is enough to make the difference between whether the star is sub- or super- 
Eddington at around 100M Q . This measure has been invoked in the spherical case as an important factor in determining 
the final mass of the star (Omukai & Palla 2003). While we regard the question of mass- limits due to feedback as being 
more complicated than a comparison of the protostar's luminosity to the Eddington value (Tan & McKee 2003; Paper 
II), Figure 6b shows that if such criteria are to be used, then the dependence on the geometry of the accretion must be 
accounted for. 

We calculate the spectra, and in particular the H-ionizing photon luminosity, S, of the various components as follows. 
When the accretion flow is optically thin, the luminosity from the stellar surface is AnrlaT^ l7 where T c ff, l is given in 
equation (B6). When the accretion flow is opaque, the effective temperature and radius of the photosphere are determined 
as described in §B.2. The assumption of a blackbody spectrum for the stellar spectrum is quite accurate for calculation of 
fluxes of hydrogen-ionizing photons (Schaerer 2002), though not for photons that ionize He. Fortunately it is the feedback 
effects from H-ionization that are much more important (Paper II). The contribution from the boundary layer accretion 
luminosity is calculated assuming the radiation emerges with a blackbody spectrum from the upper and lower surfaces 
of an annulus that extends radially from the stellar surface by a distance equal to the height of the disk's photosphere 
(estimated to be about 1.5 density scaleheights, h, from the disk model of §5; a more accurate treatment is given in Paper 
II). Finally the contribution from the inner accretion disk (r < 10r*) is included, using the surface temperatures predicted 
by the radial disk model (§5). 

The ionizing luminosities for different models are shown in Figure 7. The contribution from the protostellar photosphere 
is usually dominant. The differences between the low and high rotation parameter cases are very large, reflecting the 
different photospheric temperatures. These differences will prove crucial for determining mass limits to the star formation 
process from protostellar feedback (Paper II). Note that in the fiducial model the ionizing luminosity shows a very rapid 
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Fig. 5. — Protostcllar (solid lines) and photosphcric (dotted lines) radii as a function of protostellar mass. The dashed line in each plot is 
the zero age main sequence radius (Schaerer 2002). (a) Top left: Spherical isothermal accretion at a rate 4.4 X 10 — 3 Mq yr — 1 corresponding 
to T = 1000 K. The squares are the results of Stahler et al. (1986) and Palla & Omukai (2001). (b) Top right: Fiducial spherical isentropic 
accretion model, (c) Bottom left: isentropic case with small amount of rotation: /k C p = 0.05. The abrupt decline of the photosphere between 
30 and 40 Mq reflects the development of an inner photosphere as the density declines (see Figure 3). (d) Bottom right: Fiducial isentropic 
case with /kep = 0.5. 



increase from m» ~ 20Mq to to* ~ 40Mq, reflecting the dramatic contraction of the star and the rapidly increasing 
internal luminosity. For these accretion rates (K' = 1) the ionizing luminosity cannot increase much more quickly than 
this model, which applies to any rotating model in which the stellar photosphere coincides with its accretion surface. 
Thus without having completed the more detailed calculations of the effects of feedback on the accretion flow (Paper II) , 
we expect that feedback processes will not become important until to* is at least ~ 3OM , and this is an anticipated 
lower-limit to the mass of the first stars. 

8. CONCLUSIONS 

Recent numerical studies have indicated that the initial conditions for primordial star formation are dense, massive 
gas cores in approximate hydrostatic and virial equilibrium. These physical properties are set by the microphysics of H 2 
cooling and not by the initial cosmological density perturbations. We have presented a theoretical model for star formation 
from these cores, basing our fiducial case on the core formed at the end of the simulation of Abel, Bryan, & Norman 
(2002). Following the philosophy of Stahler, Shu, & Taam's (1980) study of low-mass, contemporary star formation, our 
approach has been to identify and isolate the various key stages that can be treated analytically and semi-analytically, 
and then combine them into a coherent picture of the whole formation process. This will be extended to include feedback 
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Fig. 6. — Top panel: Evolution of bolomctric luminosities with protostcllar mass for the fiducial iscntropic model with rotation /Kep = 0.5, 
K' = 1, and a BB = 0.01. Total radiated bolomctric luminosity (solid line) including contribution from inner accretion disk, stellar accretion 
luminosity including from direct spherical accretion and boundary layer accretion (dashed line), internal protostellar luminosity (long-dashed 
line), and accretion disk luminosity from r < 10r„ (dotted line) arc shown. The Eddington luminosity is shown by the dot- long-dashed line. 
The combination of a declining accretion rate, increasing stellar radius, and increasing stellar mass at later stages when the protostar is 
accreting on the main sequence, lead to approximately constant accretion luminosities. Note the total luminosity remains sub- Eddington over 
the entire evolution (m* < IOOOM0 . Bottom panel: Evolution of total bolometric luminosities for fiducial iscntropic models with K' = 1 and 
etas = 0.01, and with rotation /x e p = 0,0.05,0.5 (dashed, dotted, solid lines, respectively). Note that the spherical accretion case (/Kep = 0) 
exceeds the Eddington limit at m* ~ IOOMq, but the more realistic disk accretion models are slightly sub- Eddington. 



processes in Paper II. 

We have described the rate of collapse of the cores as a function of the entropy parameter of the gas, K = P/p 1 , and 
the amount of mass that has already collapsed. The mass of gas in hydrostatic equilibrium inside the point at which 
T <~ 300 K, the minimum temperature due to H 2 cooling and n ~ 10 4 cm~ 3 , the critical density of H 2 rotational-vibrational 
line cooling, yields total core masses of about 1O 3 M . These cores collapse from the inside-out to form a protostar, which 
then grows rapidly in mass, with an accretion rate m* ~ 0.017(m*/MQ) _3/ ' 7 Mq yr _1 . 

We have developed a simplified method for modeling protostellar evolution and applied the appropriate accretion rate 
for primordial protostars. The method allows for the treatment of accretion of gas with angular momentum. Specifically 
we assume angular momentum is conserved inside the sonic point of the flow. Using a realistic degree of rotation for the 
initial gas core most of the accretion occurs via a disk and conditions at the protostar rapidly become optically thin, in 
contrast to the spherical case. This means that the radiation field that the accretion envelope is exposed to is significantly 
hotter so that ionization and FUV radiation pressure feedback may become important at much smaller stellar masses 
than in the spherical case (Paper II). 
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Fig. 7. — Top: Hydrogen ionizing photon luminosities, S, as a function of protostellar mass for the fiducial iscntropic model with rotation 
/Kep = 0.5, K' = 1, and a aa = 0.01. The total ionizing luminosity (solid line) and contributions from the protostellar surface (long-dashed 
line), boundary layer (dashed line), and accretion disk from r < 10r» (dotted line) are shown. Bottom: Total hydrogen ionizing photon 
luminosities for the fiducial spherical iscntropic model (dashed), and with rotation /xop = 0.05,0.5 (dotted, solid). The zero age main 
sequence ionizing luminosity from the models of Schaerer (2002) is shown by the long-dashed line. This figure shows the crucial importance 
of the initial core angular momentum on the evolution of the ionizing luminosity: rotating systems form protostars that achieve quite high 
ionizing luminosities at relatively early stages in their evolution, compared to non-rotating, spherically-accreting systems. 



In order to determine the emission spectrum of the protostar, we evaluated the spectrum of the radiation emitted from 
the stellar surface in both the optically thick and thin cases, and we included the emission from the inner accretion disk, 
including the boundary layer. In modeling the disk, we allowed for the dissociation/ionization energy of the accreting 
gas, which can significantly lower the temperature of the disk. These predictions will be used in Paper II to estimate the 
importance of various radiative feedback processes. 

One may ask how the initial conditions and thus accretion rates of primordial star formation differ from those of 
contemporary star formation. There are several major differences: (1) the mass scale of cores is higher, being set by the 
microphysics of H2 cooling, rather than CO and dust cooling; (2) the gas temperatures are higher, which follows from (1) 
and virial equilibrium, so that at the onset of dynamical collapse the quasi-equilibrium state is denser, the free-fall time 
shorter, and the accretion rate larger; (3) because the accretion rate is larger, primordial protostars are larger and reach 
the main sequence at higher masses; (4) correspondingly, support by deuterium burning is not particularly important, in 
contrast to protostars today; (5) thermal pressure is dominant in the initial gas cores, which is not true for contemporary 
massive cores (McKee & Tan 2002), primarily because metal-rich, contemporary cores are permeated by dynamically 
important magnetic fields and are able to radiate away most of their thermal energy, leaving themselves in a state of 
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non-thermal pressure support that is supersonically turbulent; and (6) the temperature structure is such that the gas 
temperature decreases with radius, leading to accretion rates that decline with time, again opposite to the evolution of 
contemporary massive cores. 

The contraction to the main sequence (after the star is approximately older than its Kelvin-Hclmholz time), marks a 
major transition when radiative feedback from the protostar becomes much more important. Thus, even before a detailed 
analysis of these processes (Paper II), we anticipate that the masses of the first generation of stars must have been at 
least as great as <~ 30M Q . 

We thank T. Abel, V. Bromm, R. Ccn, B. Draine, J. Goodman, A. Locb, E. Komatsu, D. McLaughlin, K. Omukai, 
and J. Ostriker for helpful discussions. JCT is supported by a Spitzer-Cotsen fellowship from Princeton University and 
NASA grant NAG5-10811. The research of CFM is supported by NSF grant AST-0098365 and by a NASA grant funding 
the Center for Star Formation Studies. 



APPENDIX 

A. ENERGY ACCRETION ONTO A PROTOSTAR 

Let 

e= -v 2 + e t h + ej (Al) 

be the total energy per gram of gas, where e t h = \kJ '/ fi is the thermal energy per unit mass and ej is the internal energy 
per unit mass, including chemical binding energy. The zero point of e/ is determined by the boundary conditions. To 
relate this to the work of Nakano et al. (1995, 2000), we write ej = "fj/mn- If the gas is fully molecular, then since it 
takes 4.48 eV to dissociate H 2 , 13.6 eV to ionize H, and 79 eV to fully ionize He, the maximum value of is = 16.8 
eV for a primordial gas with 0.079 He per H. 

Conservation of energy, including radiation but neglecting viscosity and magnetic fields, is given by (see eqs. 94.15b 
and 96.17 in Mihalas & Weibel-Mihalas 1999) 



— (pe + u rad ) + V ■ 











+ F 


= Pg ■ v + pe nuc , 









(A2) 



where u ra( j is the energy density of radiation, P g is the gas pressure, F is the flux of radiation measured in the frame of 
the protostar, and pe nuc is the rate of nuclear energy generation per unit volume. This equation is accurate to order v/c. 
Now, the rate of work by the gravitational field can be written 

pg-v = -V-{pv<P)-—(p<t>)+p-^ (A3) 



(Shu 1992), so that equation (A2) becomes 

d 

— [p(e + ((>) + u rad ] + V ■ 



l>v [ e + ^ + <j> ) - F 



= P-Qi + P^nuc- (A4) 



The rate at which the protostar accretes energy can be evaluated by integrating equation (A4) over the volume of the 
protostar. Let E g be the total gas energy, etc. Then, since the gravitational energy satisfies 

dE B 



dt 

(Shu 1992), the total energy E = E g + E gr + E r obeys 

+ J dA-pv(e +L = E nuc , (A6) 

This equation is valid for an arbitrary mass distribution. The volume of integration must include all the self-gravitating 
mass, so it is valid only for r > r*. 

We now assume that the protostar is approximately spherical, so that for r > r* the potential is 

, Gm* 1 2 

<t> = — = ~2 V ^ (A7) 

where vg is the free-fall velocity. We relate the velocity of the gas to the free-fall velocity by 

v 2 

fk = -2 • (A8) 
v s 

Inside the accretion shock, fk — since we assume that the star is not rapidly rotating; in the boundary layer of an 
accretion disk, ~ 1/2; and in free-fall collapse, fk = 1. Denote the average over the surface of some quantity x by (x), 
so that if the protostar accretes matter directly via free-fall collapse at a rate direct and via disk accretion at a rate 
m^disk, then 

(*^) — ~ (^7-*, direct^direct ~i~ disk^disk) • (A9) 
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Under the assumption that the ratio of specific heats for the accreting gas is 5/3, equation (A6) for the energy becomes 



dE 
~dt 



-eth + e/--(l-/ fe )4)-L 



to* (w) — L + E n 



(A10) 



where to* is positive for accretion and w is the enthalpy plus the potential energy per unit mass. Note that in contrast 
to Nakano et al. (1995, 2000), we do not include nuclear binding energy in our expression for the energy, and as a result 
the term E nuc appears on the right-hand side 1 of equation (A10). 

Since the energy of the protostar is dominated by its interior, we are free to choose any surface layer at which to evaluate 
dE/dt. Equivalently, we can evaluate equation (A4) for steady flow in a region in which e nuc = 0. Either way, we conclude 
that the luminosity is given by 

L = L Q + to* ^|eth + ej - ^(1 - fk)v^j = L Q + m*{w), (All) 

in terms of Lq, the emergent luminosity. This equation is similar to equation (96.18) of Mihalas & Weibel-Mihalas (1999), 
except that (1) we are focusing on the surface layers so that there is no nuclear energy generation and (2) we have not 
assumed spherically symmetric flow. 

A.l. Energy Flow in an Accretion Disk 

It is customary to model accretion disks with an anomalous viscosity r], although in many cases the actual transport of 
energy and angular momentum is mediated by magnetic fields (Balbus & Hawley 1998). The energy flux due to viscous 
transport is — v • cr', where a' is the viscous stress tensor. The energy conservation equation (eq. A4) becomes 

d<j> 

pv [ e + — + 1 



«rad 



P^nuc 



(A12) 



(see Landau & Lifshitz 1959). In a steady, axisymmetric, thin (so that v\ = Gm/i 



and V • F ~ dF/dz) disk in 



which v 2 -C v 2 , (so that v 2 ~ v\ and a', is the dominant stress) and in which there are no nuclear reactions, this equation 



reduces to 



ld_ 

r dr 



pv r 
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dF 
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(A13) 



where again we have assumed that the ratio of specific heats of the gas is 5/3. Integrating through the thickness of the 
disk, we have 



2F = 
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r dr 



TO 

27 



,-eth + £/ - 



1 Gm \ 
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where to = — 2nr pv r dz is the accretion rate (defined to be positive for accretion) and W r , 
integrated stress. We have defined 



i r 



pe dz 



(AW) 

= I X oo a U dz iS the 



(A15) 



as the vertically averaged energy per unit mass through the disk, where S c = J °° pdz is half of the surface density of the 
disk. 

From the cf) component of the equation of motion, one finds that for a Keplerian accretion disk 



W r4> = - 



rhO, 
2^ 



1 



1/2' 



(A16) 



(Shakura & Sunyaev 1973), where Q is the angular velocity and ro is the radius at which the stress (<r^ = r/rdfl/dr, 
where rj is the viscosity coefficient) vanishes. Tan & Blackman (2004) have argued that a disk dynamo will create a large 
enough magnetic field that the magnetorotational instability (MM) can be effective in transporting angular momentum 
in the disk. If, however, the magnetic field in the disk is too weak for the MM to occur, two caveats must be kept in mind. 
First, in the absence of the MM, it is possible that the dominant mode of angular momentum transport will be due to 
gravitational torques in the disk, and it is not clear that the effect of these torques can be described even approximately 
by a viscous accretion disk. Second, equation (A16) is based on the assumption that the stress vanishes near the surface 
of the protostar, which is appropriate if it is slowly rotating and the transition from disk to star occurs in a thin boundary 
layer. However, for a protostar that forms through disk accretion in the absence of a magnetic field, it is not clear that 
there is a point at which dfl/dr = 0; instead, the rotation may approach that of a rigid body at small radii, becoming 
very sub-Keplerian without ever passing through a maximum. In that case, the value of ro is uncertain, but is presumably 
significantly less than the stellar radius. 

Inserting equation (A16) into equation (A14), we then find that the flux emitted from the disk is 



to d /5 _ 
inr dr \ 3 th 1 



+ ■ 



3Gtoto 
87rr 3 



1 



1/2 



(A17) 



The first term is often neglected in discussions of accretion disks (e.g., Shakura & Sunyaev 1973; Frank et al. 1995), but 
it can be important in protostellar disks. For thin disks, the thermal term e t h is negligible, but the ionization term can 
be up to an order of magnitude larger than the thermal term and can be significant. 
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B. ACCRETION DIRECTLY ONTO THE PROTOSTAR 

Gas that falls directly onto the surface of the protostar goes through an accretion shock, which we identify with the 
surface of the protostar. The accretion shock has a complex structure: the shock front, in which the gas temperature 
jumps to a high value; the post-shock relaxation layer, in which the gas cools by emitting radiation; and the radiative 
precursor, in which the gas upstream from the shock front absorbs and re-radiates the emission from the post-shock 
relaxation layer (McKee & Hollenbach 1980). The gas and radiation are thus far from LTE inside the accretion shock. 
Label the surface just outside the accretion shock front by "1" , and that just inside the postshock relaxation layer by 
"2". Stahlcr ct al. (1980) have shown that the Rosseland mean opacity of the gas between n and r 2 is small, so that the 
radiation pressure is the same at both points (P ra d,i = Pad, 2)- Since the radiation emitted by the postshock relaxation 
layer is approximately the same in the upstream and downstream directions, it follows that the ratio of the mean intensity 
J to P ra d for the radiation emitted in this layer is the same on both sides of this layer; because the layer is optically thin, 
the same is true for radiation emitted outside the layer and therefore Ji ~ J 2 . 

The luminosity in the outer layers of the protostar varies as 

L = Lq + m*w, (Bl) 
where Lq is the emergent luminosity, as seen far from the protostar, and 

^kT 1 

w=—+e I --(l-f k )4 (B2) 

(eq. All). Here fj, is the mean mass per particle at r and = v 2 jv\. For accretion onto the star, fki = 1 since the 
unshocked gas is in free-fall collapse and fk2 — since the shocked gas is nearly at rest. There is thus a jump in luminosity 
across the shock, 

In = L 2 + \f kl vl 5(^_^) +en _ £J2 . (B3) 

For free-fall collapse, L decreases beyond r\ as e t h and e/ decrease. 

B.l. Optically Thin Accretion 

To evaluate the conditions at the accretion shock, we must distinguish between the cases in which the accretion flow 
is optically thick or thin to photospheric photons; we assume that the flow is opaque to the energetic photons emitted 
by the shocked gas (cf. Stahler et al. 1980). Let the hot gas behind the shock front emit a flux of energetic photons F x 
upstream and the same flux downstream. The value of F x can be inferred from equation (B3), Aitr 2 ■ 2F X — L\ — L 2 . 
These fluxes are reprocessed into less energetic photons in the radiative precursor and the postshock relaxation layer, 
respectively. In each case, half the photons go upstream and half downstream. (Actually, Calvet & Gullbring 1998 have 
shown that the downstream flux from the radiative precursor can exceed the upstream flux due to line opacity, but we 
ignore this complication here.) As a result, a net flux F x of reprocessed photons enters the postshock gas at r 2 . These 
photons are absorbed and re-emitted by the shocked gas so that the net flux is zero. 

In this section, we are assuming that the accretion flow is transparent. We further assume that the radiation emitted 
by the hot gas in the postshock relaxation layer and by the gas in the radiative precursor is significantly more energetic 
than that emitted by the photosphere, and that correspondingly the opacity for the shock photons is significantly greater 
than the photospheric opacity (this approximation is marginal for the reprocessed photons emitted by the radiative 
precursor.) First consider the case in which the flux from the interior is negligible (Pint "C F x ). The gas at r 2 must then 
be able to radiate the flux F x into space, so that aT 2 = F x . In this case the gas inside r 2 is isothermal, so this is also the 
effective temperature T e ff, 2 . Including the effect of p nt increases the effective temperature of the gas behind the postshock 
relaxation layer to 

aT^ 2 = F int + F x . (B4) 

In the Eddington approximation the contribution of the interior flux to crT 4 is reduced by a factor 2 at the surface, since 
the radiation then occupies only half the available solid angle. As a result, we have 

o-T£ = ±F int + F x . (B5) 

Stahler et al. (1980) evaluated the temperature inside the postshock relaxation layer, at a point at which the X-rays had 
been emitted but not yet reprocessed, and they obtained a factor | in front of F x in this equation (corrected by Stahler 
1988). The calculations of Calvet & Gullbring (1998 — see below) show that the temperature structure of the shocked gas 
is actually more complicated than these simple analytic models imply, however. The effective temperature at r x includes 
the contribution of the reprocessed flux F x emitted upstream, so that 

aT c 4 ffil =P int + 2P x . (B6) 

Calvet & Gullbring (1998) carried out both analytic and numerical models of the accretion shock. In their analytic 
calculation, they divided the postshock relaxation layer into two parts, the region in which the energetic photons arc 
emitted and that in which they are absorbed, and they included the latter region in their calculation. With a flux Fi 
incident on the absorbing gas, they found 



„ 3/3 
2+-+ (q- - 



(B7) 
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where q is the ratio of the opacity of the energetic photons to the photosphcric photons. The first term is the standard 
Eddington approximation for the temperature structure of a plane parallel atmosphere carrying a flux i*i nt . If q >• 1, as 
we have assumed, then the contribution of the shock radiation to the effective temperature is 

-F t (2 + qe-vr) . 

At t = (which is inside the shock in our terminology), T 4 can be much larger than F{. Only the reprocessed photons 
penetrate into the region r < r 2 . Since an approximately equal flux of reprocessed photons enters this region from the 
radiative precursor, the solution in this region is equivalent to having Fj = 2F X . Since the energetic photons have been 
absorbed at r 2 (qr S> 1) but not photospheric photons (r <C 1), it follows that crT 2 = §(-Fint + Fi) = 5 -Pint + Fx, m 
agreement with our result in equation (B5). 

Calvet & Gullbring's (1998) solution shows that the shock is thin only if q ^> 1, as we have assumed. In this case, 
the gas emitting the reprocessed photons is hotter than the photosphere by a factor ~ q 1 ^. As a result, the radiation 
outside the photosphere does not have a blackbody spectrum; correspondingly, the analytic model for the temperature 
distribution is approximate. However, their numerical calculations show that equation (B4) for the effective temperature 
behind the shock is accurate to within a few percent. The gas between the photosphere and the shock is generally warmer 
than that at the photosphere due to the heating by reprocessed photons. 

How does the heating and ionization of the gas ahead of the shock affect the emitted luminosity? The luminosity that 
emerges from the protostar can be inferred from equation (All) 

L = L 2 -m» ( + e/ 2 - \v% ) , (B8) 



2^2 2 

where we have assumed that Lq is measured at a large enough distance from the star that To and eo are both negligible. 
To approximately allow for the non-blackbody nature of the spectrum of the radiation near the shock, we evaluate the 
ionization at the effective temperature behind the shock, e I2 = e/(T off , 2 ). 

B.2. Optically Thick Accretion 

When the accretion flow is opaque, we determine the temperature just behind the accretion shock, T 2 , from the radiation 
momentum equation (Mihalas & Weibel-Mihalas 1999, eqs. 97.2 and 97.3), 

V-P; ad = - — , (B9) 

where PJ. ad is the radiation pressure tensor and F' the flux, both measured in the comoving frame, and where k is the 
Rosseland mean opacity. The comoving flux is related to that in the protostar frame by F' = F — v(u' Tad + P r ' ad ) (Mihalas 
& Weibel-Mihalas 1999, eq. 91.17). Hence, for isotropic radiation that is in LTE (P r ' ad = a^T 4 /3) we have 

Note that v < for accretion flow, so the advection term kvT/c increases the magnitude of the temperature gradient. 
To simplify the integration, we approximate the luminosity as a constant, L ~ \{L\ + L p ), where L p is the photospheric 
luminosity. The flux is then 

F~^±^. (Bll) 

To solve for T 2 , we first guess a value of L p , which determines the photospheric temperature T p in terms of the unknown 
radius of the photosphere, r p through T p = (L p j ' Aitor 2 ) 1 ^ . We solve for the temperature structure and thus the total 
optical depth of the radiative precursor from the photosphere to the accretion shock, using the opacities of Rogers & 
Iglesias (1992) and Iglesias & Rogers (1996) for T > 7000 K and those of Lenzuni, Chcrnoff, & Salpeter (1991) for 
T < 7000 K. This is carried out at an angle of 7r/3 from the pole, so that the optical depth is a reasonable approximation 
for the mean value from the stellar surface. According to equation (Bl), L p and L 2 are related by 

L p = L 2 + ■^rhtVg^ 

5 . (kT 2 kT p \ . 

- m*(ej 2 - e Ip ). (B12) 

2 V M P / 

This value of L p is used to determine T p and the cycle is repeated until convergence is achieved. 
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